# Prepared by Tim Waring
# for publication
# 2020.07.13


# truncated normal distribution library
library(truncnorm)
library(tidyr)
library(broom)
library(dplyr)
library(ggplot2)

# sample df
hipay	<- round(rtruncnorm(30, a=0, b=275, mean = 133, sd = 57))
lopay <- round(rtruncnorm(30, a=0, b=245, mean = 103, sd = 57))
round <- c(1:30)
deck <- data.frame(round,hipay,lopay)

# test
t.test(deck$hipay,deck$lopay)

# basic approach
deck %>%
	filter(round <= 20) %>%
	do(tidy(t.test(x=.$hipay, y=.$lopay))) %>%
	pull(p.value)

# basic function
calc_cert <- function(maxround){
	deck %>%
		filter(round <= maxround) %>%
		do(tidy(t.test(x=.$hipay, y=.$lopay))) %>%
		pull(p.value)
}

calc_cert(20)

# use these settings for final run
do(deck, tidy(t.test(x=.$hipay, y=.$lopay, alternative = "two.sided", paired = TRUE, var.equal = TRUE, conf.level = 0.95)))$statistic


# vector output function
calc_certs <- function(df){

	certs <- NULL
	maxrnd <- max(df$round)
	round <- c(1:maxrnd)
	cert <- rep(1,maxrnd)
	certs <- as.data.frame(cbind(round,cert))

	for (this_round in df$round) {
		if (this_round < 2 ){
			cert <- 1
		}
		else {
			cert <- df %>%
				filter(round <= this_round) %>%
				do(tidy(t.test(x=.$hipay, y=.$lopay))) %>%
				pull(p.value)
		}
		certs$cert[this_round] <- cert
	}
	certs
}

calc_certs(deck)



# multi column output function
calc_certs_reps <- function(indf,reps){

	outdf <- NULL

	for (rep in 1:reps) {
		certs<-calc_certs(indf)
		outdf <- as.data.frame(rbind(outdf,certs))
	}
	outdf
}

calc_certs_mc2(deck,3)



# multi column output function
calc_certs_reps_in <- function(reps){

	outdf <- NULL

	for (rep in 1:reps) {
		hipay	<- round(rtruncnorm(30, a=0, b=275, mean = 133, sd = 57))
		lopay <- round(rtruncnorm(30, a=0, b=245, mean = 103, sd = 57))
		round <- c(1:30)
		indf <- data.frame(round,hipay,lopay)
		certs<-calc_certs(indf)
		outdf <- as.data.frame(rbind(outdf,certs))
	}
	outdf
}


calc_certs_reps_in(100) %>%
	ggplot (aes(x=round, y=1-cert)) + geom_smooth() + ylab("mean certainty")
ggsave("average certainty.pdf", width = 5, height = 3)


